In [1]:
import pandas as pd
#東京電力のオープンデータを読み込む
elec = pd.read_csv(
    "path to data",
#日付と時間を結合したセルに変更する
    parse_dates={'datetime': ['DATE', 'TIME']})
#使用する文字列をアルファベットに変更する
columns={
    "東京エリア需要":"tky_dmnd",
    "火力":"thp",
    "水力":"hyd",
    "太陽光発電実績":"slr"
}
#再度読み込みを行う
elec.rename(columns=columns, inplace=True)
In [2]:
elec.head()
Out[2]:
datetime tky_dmnd 原子力 thp hyd 地熱 バイオマス slr 太陽光出力制御量 風力発電実績 風力出力制御量 揚水 連系線 合計
0 2016-04-01 00:00:00 2,555 0 2,258 92 0 2 0 0 2 0 0 201 2,555
1 2016-04-01 01:00:00 2,433 0 2,151 92 0 2 0 0 2 0 0 186 2,433
2 2016-04-01 02:00:00 2,393 0 2,117 92 0 2 0 0 2 0 0 180 2,393
3 2016-04-01 03:00:00 2,375 0 2,102 93 0 2 0 0 1 0 0 178 2,375
4 2016-04-01 04:00:00 2,390 0 2,126 93 0 2 0 0 1 0 0 169 2,390
In [3]:
import plotly
plotly.offline.init_notebook_mode(connected=False)
#ラベルオプションの設定を行う
data = [plotly.graph_objs.Scatter(x=elec["datetime"], y=elec["tky_dmnd"], name="東京エリア需要")]
layout = plotly.graph_objs.Layout(
    title="東京エリア需要(1時間毎の時系列)",
    legend={"x":0.8, "y":0.1},
    xaxis={"title":"datetime"},
    yaxis={"title":"東京エリア需要(万kw)"}
)
fig = plotly.graph_objs.Figure(data=data, layout=layout)
plotly.offline.iplot(fig, show_link=False)
In [4]:
thp = plotly.graph_objs.Scatter(x=elec["datetime"], y=elec["thp"], mode = "lines", name="火力")
hyd = plotly.graph_objs.Scatter(x=elec["datetime"], y=elec["hyd"], mode = "lines", name="水力")
slr = plotly.graph_objs.Scatter(x=elec["datetime"], y=elec["slr"], mode = "lines", name="太陽光")
#火力・水力・太陽光のデータの結合を行う
el = [thp, hyd, slr]

layout = plotly.graph_objs.Layout(
    title="発電量割合(1時間毎の時系列)",
    legend={"x":1.0, "y":1.0},
    xaxis={"title":"date_time"},
    yaxis={"title":"発電量(万kw)"}
)
fig = plotly.graph_objs.Figure(data=el, layout=layout)
plotly.offline.iplot(fig, show_link=False)
In [5]:
#気象庁のデータの読み込みを行う
wh = pd.read_csv('path to data', encoding="sjis",parse_dates={'datetime': ['date', 'time']})
wh.head()
Out[5]:
datetime temp rain solor_rad wind_speed wind_direction sunshine_hours hPa humidity
0 2016-04-01 01:00:00 13.3 0.0 0.0 1.4 南西 0.0 1014.8 76.0
1 2016-04-01 02:00:00 13.8 0.0 0.0 2.6 北西 0.0 1014.6 72.0
2 2016-04-01 03:00:00 13.0 0.0 0.0 3.1 北西 0.0 1014.7 63.0
3 2016-04-01 04:00:00 12.2 0.0 0.0 2.4 北北西 0.0 1014.9 67.0
4 2016-04-01 05:00:00 11.2 0.0 0.0 1.5 北北西 0.0 1015.4 69.0
In [6]:
#気温と電力需要の相関を作成するのに必要なカラムを抽出する。結合条件は日時
dat = pd.merge(elec, wh, how="inner")
dat = dat.ix[:,["datetime","tky_dmnd","temp"]]
dat.head()
Out[6]:
datetime tky_dmnd temp
0 2016-04-01 01:00:00 2,433 13.3
1 2016-04-01 02:00:00 2,393 13.8
2 2016-04-01 03:00:00 2,375 13.0
3 2016-04-01 04:00:00 2,390 12.2
4 2016-04-01 05:00:00 2,467 11.2
In [7]:
cor_dmnd  = [plotly.graph_objs.Scatter(x=dat["temp"], y=dat["tky_dmnd"], mode = "markers")]
layout = plotly.graph_objs.Layout(
    title="電力需要と気温の関係",
    legend={"x":0.8, "y":0.1},
    xaxis={"title":"気温(℃)"},
    yaxis={"title":"電力需要(万kw)"}
)
fig = plotly.graph_objs.Figure(data=cor_dmnd , layout=layout)
plotly.offline.iplot(fig, show_link=False)
In [8]:
import plotly as py
import plotly.graph_objs as go
plotly.offline.init_notebook_mode(connected=False)

trace1 = go.Scatter3d(x=dat["datetime"], 
                      y=dat["temp"],
                      z=dat["tky_dmnd"],
    mode='markers',
    marker=dict(
        size=2,
        colorscale='Viridis', 
        opacity=0.8
    )
)

data = [trace1]
layout = plotly.graph_objs.Layout(
    title="気温と電力需要の時系列3D"
)
fig = go.Figure(data=data, layout=layout)
py.offline.iplot(fig, show_link=False)
In [9]:
#東京電力のオープンデータを読み込む
elec = pd.read_csv(
    "path to data",
#indexをdatetimeにする
    index_col="datetime",
#日付と時間を結合したセルに変更する
    parse_dates={'datetime': ['DATE', 'TIME']})
#使用する文字列をアルファベットに変更する
columns={
    "実績(万kW)":"tky_cons"
}
#再度読み込みを行う
elec.rename(columns=columns, inplace=True)
elec_dat = elec.ix[:,["tky_cons"]]
elec_dat.head()
Out[9]:
tky_cons
datetime
2016-04-01 00:00:00 2555
2016-04-01 01:00:00 2433
2016-04-01 02:00:00 2393
2016-04-01 03:00:00 2375
2016-04-01 04:00:00 2390
In [10]:
wh = pd.read_csv('path to data', index_col="datetime", encoding="sjis",parse_dates={'datetime': ['date', 'time']})
wh.head()
Out[10]:
temp rain solor_rad wind_speed wind_direction sunshine_hours hPa humidity
datetime
2016-04-01 01:00:00 13.3 0.0 0.0 1.4 南西 0.0 1014.8 76.0
2016-04-01 02:00:00 13.8 0.0 0.0 2.6 北西 0.0 1014.6 72.0
2016-04-01 03:00:00 13.0 0.0 0.0 3.1 北西 0.0 1014.7 63.0
2016-04-01 04:00:00 12.2 0.0 0.0 2.4 北北西 0.0 1014.9 67.0
2016-04-01 05:00:00 11.2 0.0 0.0 1.5 北北西 0.0 1015.4 69.0
In [11]:
#東京の電力消費データと気温のデータの結合を行う
tokyo = elec_dat.join(wh["temp"]).dropna().as_matrix()
tokyo
Out[11]:
array([[ 2433. ,    13.3],
       [ 2393. ,    13.8],
       [ 2375. ,    13. ],
       ..., 
       [ 2994. ,     5.4],
       [ 2946. ,     5.5],
       [ 2783. ,     5.4]])
In [14]:
tokyo_elec = tokyo[:, 0:1]
tokyo_wthr = tokyo[:, 1:]
In [15]:
# SVMを用いるために必要なライブラリーの呼び出し
import sklearn.cross_validation
import sklearn.svm

data_count = len(tokyo_elec)
In [16]:
# KFoldの呼び出しと定義
kf = sklearn.cross_validation.KFold(data_count, n_folds=5)

# 交差検定実施
for train, test in kf:
    x_train = tokyo_wthr[train]
    x_test = tokyo_wthr[test]
    y_train = tokyo_elec[train]
    y_test = tokyo_elec[test]

    # -- SVR --
    model = sklearn.svm.SVR()
    y_train = y_train.flatten()
    y_test = y_test.flatten()
    model.fit(x_train, y_train)
    print ("SVR: Training Score = %f, Testing(Validate) Score = %f" %
           (model.score(x_train, y_train), model.score(x_test, y_test)))
SVR: Training Score = 0.141059, Testing(Validate) Score = -0.470640
SVR: Training Score = 0.117703, Testing(Validate) Score = 0.121252
SVR: Training Score = 0.073607, Testing(Validate) Score = -0.129368
SVR: Training Score = 0.134667, Testing(Validate) Score = 0.024817
SVR: Training Score = 0.166388, Testing(Validate) Score = -0.293765
In [18]:
#予測モデルの作成
model = sklearn.svm.SVR()
y_train = y_train.flatten()
y_test = y_test.flatten()

model.fit(x_train, y_train)
Out[18]:
SVR(C=1.0, cache_size=200, coef0=0.0, degree=3, epsilon=0.1, gamma='auto',
  kernel='rbf', max_iter=-1, shrinking=True, tol=0.001, verbose=False)
In [19]:
#最小値、最大値は、それぞれ気温の最高・最低を参照する。刻み込み幅は0.01
import numpy as np
px = np.arange(tokyo_wthr.min(), tokyo_wthr.max(), 0.01)[:, np.newaxis]
py = model.predict(px)
In [20]:
ppy = pd.DataFrame(px, py)
ppy.rename(columns = {0:'temp'}, inplace=True)
In [21]:
#気温と電力需要のグラフにSVMで作成された予測モデル
pred  = [plotly.graph_objs.Scatter(x=ppy["temp"], y=ppy.index, mode = "lines")]

layout = plotly.graph_objs.Layout(
    title="電力需要と気温の関係(SVMを用いた予測モデル)",
    legend={"x":0.8, "y":0.1},
    xaxis={"title":"気温(℃)"},
    yaxis={"title":"電力需要(万kw)"}
)
fig = plotly.graph_objs.Figure(data=pred, layout=layout)
plotly.offline.iplot(fig, show_link=False)
In [22]:
#気温と電力需要のグラフにSVMで作成された予測モデルを重ね合わせる

cor_dmnd = plotly.graph_objs.Scatter(x=dat["temp"], y=dat["tky_dmnd"], mode = "markers", name="実測値散布")
pred = plotly.graph_objs.Scatter(x=ppy["temp"], y=ppy.index, mode = "lines", name="需要予測")

el = [pred, cor_dmnd]
layout = plotly.graph_objs.Layout(
    title="SVMを用いた需要予測モデル",
    legend={"x":1.0, "y":1.0},
    xaxis={"title":"気温(℃)"},
    yaxis={"title":"電力需要(万kw)"}
)
fig = plotly.graph_objs.Figure(data=el, layout=layout)
#plotly.offline.iplot(fig, show_link=False)
In [23]:
#当てはまりがわるいので、カーネル密度推計を行い、分布を立体的に確認する
from sklearn.neighbors.kde import KernelDensity
import numpy as np

dmnd_cor = tokyo[:, 0:]
#dmnd_cor
kde = KernelDensity(kernel='gaussian', bandwidth=0.2).fit(dmnd_cor)
kde_score = kde.score_samples(dmnd_cor)
kde_score
Out[23]:
array([-6.78122567, -7.41382117, -7.28689816, ..., -7.41382617,
       -6.45580608, -7.41382617])
In [24]:
import plotly as py
import plotly.graph_objs as go
plotly.offline.init_notebook_mode(connected=False)

#3D用の呼び出し、x=需要、y=気温、z=密度のグラフを書く
trace1 = go.Scatter3d(x=dat["tky_dmnd"], 
                      y=dat["temp"],
                      z=kde_score,
    mode='markers',
    marker=dict(
         size=2,
        colorscale='Viridis', 
        opacity=0.8
    )
)

data = [trace1]
layout = plotly.graph_objs.Layout(
    title="気温と電力需要の確率密度"
)
fig = go.Figure(data=data, layout=layout)
py.offline.iplot(fig, show_link=False)
In [37]:
#AERIALのデータを参照するsample code
import urllib.request
import json
with urllib.request.urlopen('http://api.hi-rezclimate.org/amsr2.py/prc?lat=35.689622&lng=139.761900&date=2016-04-01&range=0.5') as response:
        prc = response.read()
In [38]:
prc = json.loads(prc)
prc["values"][0:10]
Out[38]:
[{'lat': 36.1, 'lng': 139.2, 'prc': 9999.0},
 {'lat': 36.1, 'lng': 139.3, 'prc': 0.0},
 {'lat': 36.1, 'lng': 139.4, 'prc': 0.0},
 {'lat': 36.1, 'lng': 139.5, 'prc': 0.0},
 {'lat': 36.1, 'lng': 139.6, 'prc': 0.0},
 {'lat': 36.1, 'lng': 139.7, 'prc': 0.0},
 {'lat': 36.1, 'lng': 139.8, 'prc': 0.0},
 {'lat': 36.1, 'lng': 139.9, 'prc': 0.0},
 {'lat': 36.1, 'lng': 140.0, 'prc': 0.0},
 {'lat': 36.1, 'lng': 140.1, 'prc': 0.0}]
In [ ]: